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In order to clarify the "superradiant" conditions for the moving Josephson vortices to excite in- 
phase AC electromagnetic fields over all junctions, we perform large scale simulations of realistic 
dimensions for intrinsic Josephson junctions under the layer parallel magnetic field. Three clear 
step-like structures in the I-V curve are observed above a certain high field ( H > IT in the present 
simulations ), at which we find structural transitions in the moving flux- line lattice. The Josephson 
vortex flow states are accordingly classified into four regions ( region I ~ IV with increasing current 
), in each of which the power spectrum for the electric field oscillations at the sample edge are 
measured and typical snapshots for Josephson vortex configurations are displayed. Among the four 
regions, especially in the region III, an in-phase rectangular vortex lattice flow state emerges and 
the power spectrum shows remarkably sharp peak structure, i.e., superradiant state. Comparison of 
the simulation results with an eigenmode analysis for the transverse propagating Josephson plasma 
oscillations reveals that the resonances between Josephson vortex flow states and some of the eigen- 
modes are responsible for the clear flux lattice structural transitions. Furthermore, the theoretical 
analysis clarifies that the width of the superradiant state region in the I-V characteristics enlarges 
with decreasing both the superconducting and insulating layer thickness. 



PACS numbers: 74.50. +r, 74.60. Ge, 74. 80. Dm, 85.25.Cp 



I. INTRODUCTION 

Intrinsic Josephson effects (IJE's) in highly anisotropic 
high temperature superconductors (HTSC's) such as 
Bi2Sr2CaCu208 is a subject currently under intensive 
investigation Q| . This is because single crystals of these 
materials naturally form a stack of many atomic-scale 
Josephson junctions, i.e., intrinsic Josephson junctions 
(IJJ's) and therefore strong coherence between junctions 
are expected. Especially, if in-phase oscillation over all 
junctions is realized it will show superradiant properties 
such as remarkable increase of output power and decrease 
of lincwidth [2-7], and its realization will give a great im- 
pact on technological applications to optical communi- 
cations in the near future. Therefore, in this paper, we 
explore such superradiant conditions in Josephson vortex 
flow states in IJJ's. 

In the case of the single long Josephson junctions 
(SLJJ's), fluxon dynamics coupled with transverse prop- 



agating Josephson plasma oscillation have already been 
extensively investigated |8), ||, M], and SLJJ's have 
been employed as useful flux flow oscillators. In SLJJ's, 
the fluxons are aligned like one-dimensional chain along 
the junction plane and unidirectionally driven under the 
bias-current. This sliding motion of the fluxon chain ex- 
cites the transverse propagating Josepshson plasma os- 
cillations inside the junction, and when their moving 
speed matches the Josephson plasma mode velocity, I-V 
characterisitcs show a step like resonant structure, where 
the electromagnetic (EM) wave radiation power becomes 
maximal fioj, 0. 

On the other hand, IJJ's have a longitudinal propagat- 
ing Josephson plasma mode along the stacked direction 
and its dispersion relation are shown to be very different 
from that of the transverse mode, i.e., the longitudinal 
one is almost flat while the transverse one is very disper- 
sive. Thus, IJJ's can have many propagating modes with 
different dispersions as mixtures between transverse and 
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longitudinal modes, and IJJ systems with finite stack- 
ing height (employed in experiments) have the standing 
waves for those mixed modes along the stacked direc- 
tion |Q , H . Accordingly, the I- V characterisitics and the 
properties of EM radiation in IJJ's are clearly expected 
to exhibit a rich variety in contrast to SLJJ's. Actually, 
so far, the fluxon dynamics in IJJ's have been investi- 
gated in several papers[2-7], Q to understand their dy- 
namics and related EM radiation, and complex resonant 
behaviors with some special modes have been reported. 
However, the full understanding for the behavior of the 
fluxons and Josephson plasma oscillations has not been 
achieved. Especially, the clear in-phase superradiant con- 
ditions in IJJ's have not yet been established, and some 
questions still now remain concerning the stability of the 
superradiant flow state and its emergent region with re- 
spect to experimental variables as the magnetic field and 
the transport current. Therefore, in this paper, we report 
results of the first direct large scale simulations for the 
fluxon dyanamics in IJJ's using the coupled Sine-Gordon 
equation and clarify the superradiant conditions. Fur- 
thermore, we pay attention to all dynamical flux flow 
states observed in various current ranges in the I-V char- 
acteristics and reveal that the I-V characteristics show 
three clear step-like structures corresponding to flux lat- 
tice structural transitions. Thus, we distinguish flux flow 
states into four regions at these three step-like structures, 
and display the typical flux flow lattice configuration in 
each of the regions. The names region I ~ IV ( in in- 
creasing order of the current ) will be assigned to the 
four regions. In order to identify the natures of the EM 
radiations in the four regions, the electric field oscilla- 
tions at the sample edge, which characterize the proper- 
ties of the emitted electromagnetic wave, are measured 
and their typical power spectrum in each region is given. 

Here, let us briefly describe the characters of the flux 
flow states observed in the four regions. The notable 
points are remarkable differences in the flux lattice con- 
figurations and power spectra for EM field oscillations in 
those four regions. In region I, fluxons flow randomly 
without forming a clear lattice structure. This state is 
a steady state when the applied current is low enough, 
and the electric field oscillations at the edge exhibit a 
broad-band power spectrum. In region II, wavy aligned 
flux flow states appear. From this region, the vortex 
flow states begin to show resonances with some of the 
transverse propagating modes. The vortex alignments 
develop along the stacked direction, and the power spec- 
trum has peak structures at a characteristic frequency 
and its higher harmonics. In region III, a perfect in- 
phase flux flow states emerges. The vortex alignment 
becomes perfect along the c-axis and fluxons form a rect- 
angular lattice. The line- width of the power spectrum 
is very sharp and the power at the peak position is also 
very large, i.e., the superradiant state. In region IV, 
triangular like lattice flux flow states are observed. Its 
Josephson frequencies due to the flux flow voltages are 
larger than the frequencies of any plasma modes. Conse- 



quently, it is expected that the transverse plasma modes 
does not affect the formation of flux lattice structures 
and therefore energetic favorable triantular like lattice 
flow states appear. In order to understand the origin 
of the transitons between these four flux flow states, we 
perform an eigenmode analysis on the employed coupled 
Sine-Gordon equation. As a result, we are able to iden- 
tify that the resonances between three modes and the flux 
flow states leads to three clear step-like structures in the 
I-V characteristics. Thus, from this theoretical analysis, 
we can predict the width of the four regions, especially, 
the superradiant state (region III), in the I-V character- 
istics. Also, the analysis clarifies that the superradiant 
state region enlarges with increasing a parameter ratio, 
a 2 

-^p-, where s and D are, respectively, superconducting 
and insulating layer thickness and \ a b is the penetration 
depth in the ab-plane direction. This result clearly indi- 
cates that IJJ's can have remarkably wide superradinat 
region in I-V characterisitics. 

The outline of this paper is as follows. The model 
equation employed is briefly introduced in section II, and 
the numerical simulation results are shown in section III. 
Section IV is devoted to eigen mode analysis and discus- 
sions. 



II. MODEL EQUATION 

We first explain our model. We perform the simula- 
tions under a configuration in which the transport cur- 
rent flows in the direction perpendicular to the junctions 
(|| z-axis) and the magnetic field is applied along the 
junctions (|| y-axis) as shown in Fig.l. 



Transport Current 




►X 

FIG.l A schetch of the intrinsic Josephson junction. 
The magnetic field and the transport current are applied 
along the junction plane (y-axis) and the stacked direc- 
tions (z-axis) , respectively. The Josephson vortices move 
in the -x direction. 
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The Josephson vortices move parallel to the — ir-axis in 
this configuration. Assuming that the superconducting 
layers have an atomic-scale thickness, we neglect the vari- 
ation of physical quantities along the z-axis inside the 
superconducting layers. Then, the vector potential is de- 
fined only inside the insulating layers as e (x, t)(= 

jf^ 1 dzA z (x,t)), where I stands for the layer index, 
and the electromangtic fields Ef +l e (x, t) and Bf +1 e (x, t) 

are, respectively given by Ef +1 e (x,t) 



+1. 

and C u (i,t) - ^ 
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assume Maxwell's equations in the present system as fol 
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where t' — Lu p t and x' = X c x with uj p and A c being, re- 
spectively, the Josephson plasma frequency (u> p = ) 
and the penetration depth in the c-axis direction. The 
parameter /?(= 4 ^^ c ) is related the McCumber param- 
eter (3 C as (3 C = 1//3 2 . It is noted that eq.(5) gives a 
flat dispersion relation for the longitudinal mode along 
the c-axis while it leads to the dispersive one along the 
junctions. 
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where pi and e c are, respectively, the charge density at 
£-th superconducting layer and the dielectric constant of 
the insulating layers. The current, j Z £ + u in eq.(l) is 
given by a sum of Josephson and quasiparticle currents 
as j c sinPi+i/ + aE z £+i^£ , with j c and a being the critical 
current density and conductivity of the quasiparticles, re- 
spectively. Pi+i,i is the gauge invariant phase difference 
defined as Pg+i.e = ®e+i - t - ^A i+1 j with 6 e being 
the phase of the superconducting order parameter of ^-th 
layer. The time and spatial derivatives for Pe+ij yield 
the generalized relations, 
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U*i+i-3*i) + -T-By i+u , (4) 



where A a j, are the penetration depth associated with the 
current j x t in the ab-plane and /x is the Debye length ]13[ ] 
[Q. In deriving eq.(3) and (4), we employ the equation 
= + wnic h is analogous to the Lon- 

don equation and the London equation for the in-plane 



current, j x e = — 



4ir\ a 



r{^V x 6i-A x i). In the case where 



jf +1 — jf ^ under an external magnetic field the first 
term on the right hand side of eq. (3) gives only small con- 
tribution of 0{(x 2 /}? ab ) to the dynamics of Pt+i t, which 
means that the incomplete screening effect for charges in 
the superconducting layers, namely the coupling between 
junctions due to the charge neutrality breaking effect [jl3| , 
fl4| may be neglected in the presence of an external mag- 
netic field. Then, in this approximation the equation for 
P e+h e is derived from eqs.(I)-(4) as follows @, @, 0] 
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III. RESULTS OF NUMERICAL EXPERIMENTS 
AND DISCUSSIONS 



In this paper, we numerically solve eq.(5) for a finite 
size IJJ composed of 20 junctions in which the ac cross 
section ( the xy-plane in Fig.l) has 2-dimensional (2D) 
rectangular shape in the presence of an applied field (|| y- 
axis) and a transport current (|| z-axis). The transport 
current and the magnetic field which includes the self- 
field are introduced through the boundary conditions at 
the top and bottom layers and the edges of all the lay- 
ers. To perform the simulations we divide the system 
into 400 meshes along the x-direction. Each mesh is as- 
sumed to be 0.1 X a b wide along the ir-direction. Then the 
system used in the simulations has the area of 10^m x 20 
junctions if we assume A c = 125/ito and 7 = 500 for 
the anisotropic parameter, which are typical parameter 
values in Bi2Sr2CaCu20s. The value of the material pa- 
rameter (3 is chosen as j3 = 0.1(/3 C = 100), which is also 
a typical value in Bi2Sr2CaCu20g. The results shown 
in this paper are qualitatively size-independent. We con- 
firmed this by performing the simulations for system sizes 
up to 30/im x 80 junctions in real scale. 

Let us now present the results of numerical simula- 
tions. The applied magnetic field used in the following 
calculations is fixed to a value, 2T. The results for lower 
fields (< IT), which are qualitatively different, will be 
published elsewhere. Figure 2(a) shows the I — V char- 
acteristics of IJJ. As seen in this figure, three clear steps 
appear at the voltage values indicated by arrows in the 
I — V curve. To understand the origin of the steps, we 
plot the time and spatial averages of the Poynting vec- 
tor, < P_ x >=< j^g JdxJ2 e E^ +he B y e+i e >, and the 
Josephson current, < J z >=< -i f dx^2 e smPg + i^ > in 
Fig. 2 (b) and (c), where <> stands for the time average 
and S is the area of the computational region. Note that 
< P- x > represents the energy flow of the transverse 
EM field and < J z > gives a measure of the viscosity 
against the driven flux-flow. It is seen that the anomalous 
structures also appear in these averages, i.e., the step-like 
structure in < P_ x > and the peaks in < J z >, at the 
same voltage values as in the I —V characteristics. Such 
anomalies are commonly observed in our simulations in 
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the case where the magnetic field is greater than IT. 
Thus, it is found that the vortex dynamics are directly 
coupled with the transverse EM field and the vortex vis- 
cosity drastically changes at the step like structures. Es- 
pecially, it is noted that the viscosity of the Josephson 
vortices takes maximum values at the steps. From this 
fact one may expect that the collective vortex dynam- 
ics changes at the steps. In the following, we, therefore, 
separate the voltage region into four from I to IV at the 
steps as shown in Fig. 2(a) and investigate the collective 
properties of each flux-flow state separately. 




2000.0 4000.0 
V (a.u.) 
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FIG. 2 (a) The I-V characteristics, (b) The averaged 
poynting vector in x-direction < P- x > vs. V. The inset 
of (b) gives the enlarged view of the region indicated by 
the inclined arrow, (c) The averaged Josephson current 
vs. V. 
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FIG. 3, The snapshot of the vortex center positions 
along with sinPg+i^ distributions in all junction sites 
and the power spectrum of the electric field at the left 
edge in the regions from I to IV: (I) I/J c = 0.20, (II) 
I/J c — 0.52. Note that the Josephson plasma frequency 
at the zero temperature, u> p = ^ , corresponds to 
ijj = 130 in the figures for the power spectrum. 

Let us study the flux- flow state in region I first. In 
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Fig. 3 (I) we show a snapshot of the distribution of the 
Josephson current, sinP^+i^, and the centers of the 
Josephson vortices ( indicated by dots ) together with the 
power spectrum of the electric field at the left edge of the 
sample, S„ = J Q T Et+U (0, t)e l ^dt\ 2 . As seen 

in this figure, the vortices are distributed almost ran- 
domly and the power spectrum shows to~ 2 dependence 
above some frequency, which indicates the irregular mo- 
tion of vortices. We confirmed this chaotic behavior of 
moving vortices by monitoring directly the time devel- 
opment of the vortex positions too. From these results 
one may conclude that the highly viscous flux-flow state 
is realized in region I. The fluxons are easily pinned by 
any type of imhomogeneities, e.g., crystal imperfections, 
sample boundary etc. in this low driving-current regime. 
In the present system without imperfections the sam- 
ple boundaries including the top and bottom junctions 
work as effective pinning centers, which is seen from the 
fact that the number of vortices in the top and bottom 
junctions is much smaller than that in other intermedi- 
ate junctions. It is also noticed that the power spectrum 
forms a broad band in the low frequency region. Thus, 
the chaotic viscous flux motion generates a broad-band 
oscillating EM field in region I. The intensity of the emit- 
ted EM field increases with increasing voltage values in 
this region. 

Let us next turn to the region II. As seen in the snap- 
shot given in Fig. 3 (II), a wavy-chain-like vortex distri- 
bution is stabilized in this region. This result indicates 
that a correlation between vortices is developed in re- 
gion II. Accordingly some characteristic frequency and 
its higher harmonics appear inside the broad band spec- 
trum of the EM field. These frequencies are related to 
the quasi-periodicity of the moving vortex lattice. It is 
also noted that the frequency of the highest peak seen 
in the figure is about eight times larger than the Joseph- 
son plasma frequency at zero temperature, to p {= )■ 
This frequency corresponds to about 4THz in the case of 
A c = 125/xm and e c = 25. 

When the voltage is further increased and reaches the 
value at the boundary between the regions II and III, 
the moving vortices form a rectangular lattice as shown 
in Fig. 3(111). This regular distribution of the vortices is 
stable and maintained throughout the region III. It is 
seen that sharp peaks appear in the power spectrum of 
the EM field. The frequencies giving the peaks are de- 
termined by the flux flow velocity and the period of the 
lattice. Note that the phase differences of all the junc- 
tions make perfect in-phase motion in this state, that 
is, the superradiant state is realized in the region III. In 
this superradiant state one may expect the strong radi- 
ation of the coherent EM field from the junctions. The 
frequencies of the emitted EM field will be tunable in a 
wide range by the applied magnetic field or the voltage 
since the region III is broad as seen in the I — V curves. 
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FIG.3 (III) I/Jc = 1.16, (IV) I/J c = 1.90. 

The rectangular lattice of the moving vortices becomes 
unstable at the boundary between the regions III and 
IV. The moving vortex lattice appearing in the region IV 
is a modulated triangular one as depicted in Fig. 3 (IV). 
The power spectrum given in Fig. 3 (IV) indicates that 
the intensity of the EM field generated in this state is 
greatly reduced from that in the region III. This result is 
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easily understood in the following way. The oscillatory 
components of the tunneling current on two consecutive 
junctions cannot be in-phase in a triangular-like lattice. 
As a result, the EM field generated at the junction edges 
on consecutive junctions interfere destructively with each 
other and, therefore, the intensity of the power spectrum 
drastically decreases in this region. 



IV. EIGEN MODE ANALYSIS AND 
DISCUSSION 



To understand the origin of the transitions in the mov- 
ing vortex state we perform a linear analysis in the fol- 
lowing. The linearized equation obtained from eq.(5) 
in the case of (3 — yields the eigenvalues, uj n (k x ) — 

{1 + kl/[l + ( 1 - cos ^ )] } 1 ' 2 , and the eigenfunctions, 

P w (x',t') = e^'e^'f^ with fW u = g n e Mn / N 
for the system with N junctions, where n specifies the 
eigenmode, i.e., n = 1, ■ ■ ■ , N —1. These eigenmodes may 
be interpreted to be the transverse Josephson plasma 
modes in finite IJJ composed of N junctions without 
vortices. In the inset of Fig. 4(a) we plot the spa- 
tial variation of Re(/j™^ e ) for four values (n=0,l,2,and 
19) of n in the case of N = 20 , which represents 
the standing waves along the c-axis. In the presence 
of an external magnetic field, H, and a DC voltage, 
V, eq.(5) has an approximate solution corresponding 
to the phase-locked flux-flow state as Pt+\ t t(x' , t') — 
H'x' - V't' + he+i,e{x',t'), where H' = and 
V' = j^j-V- Here, the function he+i t t(x' ,t') is as- 
sumed to be an small oscillatory component. In this 
case, h£+i^i(x' , t') satisfies approximately the equation, 

[(1 - ^A< 2 ))<9 2 - <9 2 ]^ +M = - sm{H'x' - V't 1 ), where 
A< 2 ) is defined by A^h e+1 j = ^ +2 ,£ + i + /i«-i-2/i£ +M 
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This equation has solutions showing resonances 

H' /{[I + (1 - cos if)]} 1 / 2 . This resonant 
behavior originates as a result that the Josephson oscil- 
lations induced by the moving flux-line lattice match the 
transverse Josephson plasma frequencies. Note that the 
resonant frequencies coincide with the eigen-frequencies 
w n (k x ) for k x = H' ^> 1. Figure 4(a) shows the matching 
voltages (V' = ui n (k x — H')) as a function of H in the 
case of N = 20. In Fig. 4(b), we superpose the vertical 
lines indicating the matching voltages at k x = H(= 2T) 
on the I — V curve. As seen in this figure, the voltages 
matching with the modes A (n = 0) and B (n — 1) coin- 
cide with the 3rd and 2nd steps and the lowest matching 
voltage is also located close to the 1st step in the I — V 
curves. These analyses indicate that the structural tran- 
sitions in the moving flux-line lattice are caused by the 
matching of the phase oscillations arising from the vortex 
motion with that of the eigen-modes of the system. 
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FIG. 4, (a)The magnetic field dependence of the match- 
ing voltage between the flux flow states and all transverse 
eigenmodes. The inset shows the variations along c-axis 
of Re{f n \ +hi ) of A(n = 0), B(n = 1), C(n = 2) and 

D(n = 19) modes, where f {n \ +ll = g n e Mn/N . (b) The 
I-V characterisitcs same as Fig. 1(a). The vertical lines 
indicate the matching voltage of all eigenmodes, and A, 
B, C, and D correspond to those shown in Fig. 4(a), re- 
spectively 

Here, it is noted that not all the matching voltages 
cause clear step-like anomalies in the I — V character- 
istics. However, we have directly observed the structural 
change, for example, near the voltage value at C by mon- 
itoring the vortex configuration. Since the I — V char- 
acteristics is given by a sum of the contribution from all 
the junctions, it is understood that it is not so sensi- 
tive for a slight structural change in the moving flux- 
line lattice, that is, in a case that the current distri- 
bution is not largely altered at the transition, a clear 
anomalous structure will not appear in the I — V char- 
acteristics. We also notice that the matching with the 
modes, A and B, obtained in the above linear analy- 
ses gives the precise phase boundaries of the region III. 
This fact indicates that the above analysis is sufficiently 
valid in this region, that is, the oscillatory part /i^+i/ is 
only a small perturbation and the phase-locked in-phase 
flux-flow state, which is the moving rectangular-lattice 
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state, is very stable. The width of the region III is given 



by AV = uja(H) - uj b (H) ~ [1 



sD 



V2\ab ^/1-cos(tt/JV) 



\H. 



Since 



/sD 
V2A a 



is very small in high-XI-superconductors, the 



superradiant state is expected to be very wide in IJJ's. 

Hechtfischer et al. observed a broad-band emission 
in the flux-flow state of Bi 2 Sr 2 CaCu208 single crystals, 
whose power spectrum takes a maximum around the volt- 
age at which the step-like structure appears in the I —V 
characteristics |18 . They claim that the EM waves emit- 



ted from the junctions are generated by the Chrenkov 
radiation due to the moving vortices [^8) , |h| . From our 
simuration results, we may suggest two possibilities for 
the origin of the broad-band emission. The first one is 
the Cherenkov-type radiation as claimed by Hechtfisher 
et al. Rough estimation for the present model indicates 
that the threshhold voltage for the Cherenkov radiation 
lies in region II, in which we really observe the broad- 
band spectrum in the low-frequency region as seen in 
Fig. 3 (II). It is noted that this result is consistent with 
simulation resutls by Hechtfischer et al |is|| . In this sce- 
nario, the step-like anomaly observed in the I — V curve 
of Bi-2212 may be identified with the boundary between 
the regions II and III in our simurations. On the other 
hand, it is also possible to interprete that the step core- 
sponds to the boundary between the regions I and II, 
which leads to another possibility for the origin of the ob- 
served broad-band emmission. In this second scenario, it 
is understood that the observed EM field originates from 
the vortex state in region I. Since no transverse propa- 
gating Josephson plasma mode can be coupled with such 
the slow vortex motions in this region, one may interprete 
that the broad-band spectrum reflects chaotic EM field 
oscillations due to the non-regular motion of vortices. 
Here, note that step structure in IVC at the boundary 
between the region I and II in simulations is the most 
remarkable structure and it is similar with experimental 
results . To clarify the mechanism of the broad-band 
emission, we will need more detailed theoretical and ex- 
perimental studies. 



V. CONCLUSION 

In conclusion, we performe direct large scale numerical 
simulations on Josephson vortex flow states in IJJ's and 
find structural transitions of the moving flux-line lattice 
accompanied by the step-like structures in the I-V char- 
acteristics and changes of power spectrums of AC electric 
field at the sample edge. The four typical flux flow states 
and the corresponding power spectra are identified and 
the in-phase superradiant flux flow state are found to ex- 
ist stably in a wide region in the I-V characteristics. The 
comparison between the simulation results and the eigen- 
mode analysis for the transverse propagating Josephson 
plasma modes show that those step-like structures in 
the I-V characteristics are due to resonances between 



flux flow motions and transverse propagating Josephson 
plasma modes. Furthermore, the analysis predict that 
the superradiant state region enlarges as the supercon- 
ducting and insulating layer thickness is decreased com- 
pared to the magnetic pentration depth in the ab-plane 
direction. The collective Josephson vortex flow under the 
presence of many propagating modes includes very rich 
varieties of physics, while the stable superradiant flux 
flow state will be a key feature in making successfully 
submilimeter-wave generator employing IJJ's. 
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